flASSACHUSETTS INSTITUTE OF TECHNOLOGY
A.I. LABORATORY
Artificial Intelligence October 1973
demo. No. 235
ON LIGHTNESS
Berthold K. P. Horn
The intensity at a point in ^n image is the product of the
reflectance at the corresponding object point and the intensity of
illumination at that point. Lie are able to perceive llghtneee, a
quantity closely correlated ulth reflectance. Hon then do Me eliminate
the component due to illumination from the Image on our retina? The tuo
components of image Intensity differ in their epatial distribution* A
method Is presented here which takes advantage of thle to compute
lightness from image intensity in a layered, parallel fashion.
The method is developed for a restricted class of images first used
by Land in presenting his retinex theory of color. In thle theory the
problem of color perception is reduced to one of judging black and white
lightness on three images taken in different parts of ths visual
spectrum. The method described h*re fills the need for a lightness
judging procees.
The theory has implications for potential special purpose hardware
In Image sensing devices. It should also bs of interest to cognitive
peychologiete elnce it can explain certain effects observed in the human
vleu*l system ae well ae predict net* ones. Further, the theory provides
neuro-phyelologlsts with suggsstlons about the function of certain
etructuree In the primate retina.
This report describes research done at the Artificial Intelligence
Laboratory of the Massachusetts Institute of Technology. Support
for the laboratory's artificial intelligence research is provided
in part by the Advanced Research Projects Agency of the Department
of Defense under Office of Naval Research contract N000U-70-A-
OM2-OM5.
Reproduction of thle document. In whole or in part, la permitted for any
purpose of the United States Government*
1. Revieu. 3
1.1 Theories of Color Perception. 3
1.1.1 Tri-Stlmulue Theory. *
1.1.2 Color Conversion. 4
1.2 Land's Retinex Theory. S
1.2.1 Lightness Judging. 6
1.2.2 fllni-HorJd of Mondriane. 6
1.3 Uhy Study the One-dimensional Case? 7
1.3.1 Notation. 8
1.4 One-Dimensional Method - Continoue Case. 9
1.4.1 One-Dlmenslonal Continoue Method: Details. 9
1.4.2 Normalization. 12
1.5 One-Dimensional Method - Diecrete Case. 13
1.5.1 One-Dimensional Di ecrete flethod: Details. 13
1.5.2 Selecting the Threshold. IS
1.5.3 Accuracy of the Reconstruction. 17
1.5.4 Generalizations. 17
1.5.5 Physical Models of the One-Dimensional Discrete Process* 18
2. Lightness in Two Dimensional Images. 19
2.1 Tuo-Dineneional Method - Continoue Case. 19
2.1.1 Applying the Laplaclan to a Mondrian. 19
2.1.2 Inverse of the Laplacian Operator, 29
2*1.3 Uhy one can use the Convolutions! Inverse. 23
2.1.4 Normalization. 23
2.1.5 Tuo-D i mens I ona I Continous Methodr Detal Is- ♦ - * 24
2.2 Tuo-Dimeneional Method - Discrete Caee. 25
2.2.1 Teeselatlon of the Image Plane. 25
2.2.2 Diecrete Analogue of the Laplacian. 26
2.2.3 Inverse of the Diecrete Operator. 27
2.2.4 Computational Effort and Simplification. 29
2,2*5 Tuo-Dimeneional Discrete Methods Oetatle. 38
2.2.6 Simplicity of the Inverse. 32
2.2.7 Iterative Methods of Solution, « 32
2.2.8 Convergence of Iterative and Feed-back Schemes. 33
2.2.9 Setting the Threehold. 34
2.2.18 Some Notes on This Msthod. 35
2.2*11 Dynamic Range Reduction. 36
2.2.12 A Frequency Domain Interpretation. 36
2.3 Physical Mode J e. 37
2.3.1 A Diecrete Physical Modsl. 38
2.3.2 A Feed-back Schsme for the Inverse. 48
2.4 Limitations of the Simple Scheme Presented. 42
2.5 Computer Simulation of the Discrete Method. 44
2.5.1 Form of Inveree ueed in the Computer Simulation. 4S
3. Impl icatione and Concluelone* 47
3.1 Parallel Image Proceeeing Hardware. 47
3.2 Cognitive Peychology, *7
3.3 Neuro-phyeiology. 4 *
3.4 Conclusion. **
ON LIGHTNESS 2
Acknoulegaente
I wish to thank Dr. David flarr for a numbsr of most Interesting
dlacueeione and for encouraging na to commit this msthod to papar. Ha
first observed that the primate retina appears to have Juet the right
detailed etructure to implement the functions here described.
LIGHTNESS* Definition
The relative degree to which an object reflecte light.
The Random House Dictionary
The attribute of object colors by which the object appears to reflect or
transmit more or less of the incident light*
Uabster's Seventh Nsu Collaglata Qictionary
Preview
Part 1 Is a review of the relevant Information relating to color
vleion and lightness. Thie Includes a discussion of the Land ratine*
model in a form suitable for the developments of the next part.
In part 2, Land's one-dimensional operation will be extended to
two-dimensional Images. The method depends on a layered, parallel
computation suggestive of both biological and artificial
implementations.
In part 3 some of the implications are explored and the information
on the new image processing technique is summarized.
ON LIGHTNESS 3
1. Review,
1.1 Theories of Color Perception.
Thar* has always been great interest in how we perceive colors and
numerous explanations have been forwarded (Newton 1704, Goethe) 1818,
Young 1820, Maxwell 1856, Helisholt* 1867, Hering 1875). The huaan
perceptual apparatus Is remarkably eucceeful in coping with large
variations in the i I lueinatlon. The colors us perceive are closely
correlated ulth the surface colours of the objecte viewed, deepite large
temporal and epatial differences in color and intensity of the Incident
light. Thie is eurprieing eince we cannot eenee reflectance directly.
The light inteneitg at a point in the iaage le the product of the
reflectance at the corresponding object point and the intensity of
illumination at that point - aside from a constant factor that depende
on the optical arrangment. There must then be soee difference between
theee two components of Image intensity which allows ue to dlecount the
effect of one. The two componente differ In their spatial distribution.
Incident light inteneitg will usually vary smoothly, with no
di econt inui t ies, uhlls reflectance will have eharp diecontinul ties at
edges where objecte adjoin. The rsflsctance being relatively constant
between euch edgea.
ON LIGHTNESS k
1*1.1 Trl-Stlmului Theory*
Some facto about how we see color are fairly uell established. It
appears that us have three kinds of sensors operating In bright
1 1 lull n*tl on, Hith peak aanai tivi tlee in different parts of the visible
spectrum. This is why it takes exactly three colors in additive
mixture to natch an unknoun color. I4h1 le It la a bit tricky to Measure
the sensitivity curvss of the three sensors directly, a linear transform
of these curves has been known accurately for some tine (Brlndley 19GB).
Theee curves, ceiled the standard observer curves, are sufficient to
allow one to predict color matches made by subjecte with normal colour
vieion (Hardy 193G1 .
The elmplest theory of color perception then amounte to locally
comparing the outputs of three such sensors and assigning colour on thle
basis (Young 1820, Helmholtz 18671. This however totally falls to
explain the observed color contancy. Perceived color does not depend
directly on the relative amounts of light msasured by the three sensors
(Land 1959, Lettvin 1967).
1.1.2 Color Conversion.
A number of attempts have been mads to patch up this theory under
the rubrics of "discounting of the lllumlnanf, "contraet effect
adjuetment" and "adaptation". The mors complicated theories are based
on models with large numbers of paramstsrs which era adjusted according
ON LIGHTNESS S
to empirical data IHelson 1938 & 1943, Judd 1348 & 1952, Richard* 1971).
Theee theories are at leaet partially effective in predicting human
color perception when applied to simple arrangements of stimuli similar
to those used in determining the parameters.
The parameters depend strongly on the data and alight experimental
variations will produce large fluctuations in them. This Is a phenomena
familiar to numerical analysts fitting curves to data when the number of
parameters is large. These theories are lacking in parsimony and
convincing physiological counter-parts. Lettvln has demonstrated the
hopeleesneae of trying to find flxsd transformations from locally
compared output of sensors to perceived color {Lettvln 1967).
1.2 Land*s Retinex Theory.
Another theory of color perception is embodied In Land'e retinex
model (Land 1953. 1964 & 19711. Land proposes that the three eete of
seneore are not connected locally, but instead are treated ae If they
represent points on three separate images. Processing is performed on
each euch image separately to remove the component of Intensity due to
illumination gradient. Such processing is not msrely an added frill but
ie indi sponsible to color perception in the face of the variability of
the i I lumination.
ON LIGHTNESS 6
1*2.1 Lightness Judging*
In essence a judge of lightness processes each image* Lightness is
the perceptual quantity closely correlated with surface reflectance.
Only after this process can the three images be compared to reliably
determir>e colors locally. It remains to mechanize this process.
It Mould appeal to Intuition If this process could be carried out
in a parallel fashion that doss not depend on previous knowledge of the
scene viewed* This is because colors are so immediate, and seldom
depend on one's interpretation of the scene. Colore will be seen even
when the picture makes no ssnse in terns of previous experience. Also,
color Is eeen at every point In an Image.
1*2.2 lllni-world of tlondrians.
In developing and explaining his theory Land needed to postpone
dealing ulth the full complexity of arbitrary scenes* He eelected a
particular claes of objscts as Inputs, modelled after the paintings of
the turn-of-the-century Dutch artist Pieter Cornells (londrian. These
scenes are flat areas divided into sub-regions of uniform matte color*
Problems such as those occasionsd by shadows and specular reflection are
avoided In this way. One also avoids shading; that is, the variation In
reflectance with the orientation of the surface In respect to the sensor
and the light-source. For Hondrians* lightness is considered to be a
function of reflectance*
ON LIGHTNESS 7
llondrlans ara usually made of polygonal regione mth straight aides
— for the development here however the edges may be curved. In the
world of Mondrians one finds that the reflectance has sharp
discontinuities wherever regions meet, being constant inaida each
region. The illumination, on tha other hand, varies anoothlg over the
image,
1*3 Uhy Study the One-dimensional Case?
Imagee MrB tuo-dlmeneional and usually sampled at discrete points.
For hlatoric reasons and intuitive simplicity tha results will flret be
developed in one dimension, that is with functions of one variable*
Similarly, contlnoue functlone uill be used at flret eince they aJIOM a
cleaner separation of the two componsnts of image intensity and
iiluetrate nore clearly the concepts involved.
Use Hill be Bade of analogies between the one-dimensional and two-
dimensional casss as well ae the continous and discrete ones. The final
procaee discussed for processing image intensities Is tuo-dlmenelonal
and dlecrete. A number of physical lmplenentatlons for this scheme ar9
suggested.
The proceae hIII be looked at from a number of points of view:
partial dl f ferent ial equat ione, I inear systems, four ter transforms and
convolutions, difference equations, iterative solutions, feed-back
schemes and physical models.
ON LIGHTNESS 8
1*3.1 Notation.
The following notation mi I I be used.
a* Intensity of incident illumination at a point on the object*
r' Reflectance at a point of the object*
p* Intensity at an image point* Product of •' and r'.
a t r. pi Logarithms of e\ r* and p' respectively.
d Result of applying forward or differencing operator to p.
t Result of applying threshold operator to d.
I Result of applying inverse or summing operator to t*
Simple derivative operator in one dimension*
T Contlnous threshold operator, discards finite part.
1 Simple integration operator in one dimension,
L Laplacian operator - eum of second partial derivatives.
G' Inveree of the Laplacian, convolution uith (1/2 t ) log c d/r).
D*, T*, I*, L* and Gfri Dlacrete analogues of 0, T, I, L and G.
The output It will not be called lightness since there le probably
not yet a generally acceptable definition of this term. It is however
intended to be monotonical ly related to lightness. Note that I le
ON LIGHTNESS 3
related to the logarithm of reflectance, uhile the perceptual quantity
le perhape more closely related to the square-root of reflectance.
1.4 Qne-Oimensional Plethod - Continous Case.
Land invented a simple method for separating the Image components
in one dimension* First one takes logarithms to convert the product
Into a sum* This Is followed by differentiation. The derivative uill
be the sum of the derivatives of the two components. The edges uill
produce sharp pulses of arsa proportional to the Intensltg eteps between
regione - uhi le the spatial variation of II I urn (nation uill produce only
f i nf to values overyuhere. Noh if ono diacarda all finite yatuae, one It
left ulth the pulses and hsnce the derivative of lightness* Finally one
undoee the differentiation by simple integration.
1.4.1 One-Dimensional Continous Method! Details.
Ue have the following. Let r'UI be the reflectance of the object
at the point corresponding to the image point x. Let eMx) be the
intensity at this objact point. Let p* (x) be their product, that ie,
the intensity recorded In the imags at point x. Note that e' (x) and
r* (k) %r9 poei tive,
p* <x) - e' (x) * rMx)
.«
t
1
T
dlx
T
tlxj
■ H
FIGURE 1: Processing steps in the one-dimensional continous case*
ON LIGHTNESS 11
Now let p(x) be the logarithm of p' (x) and bo oni
ptx) - six) + r(x)
Note that six) is continous and that r(x) has some finite
diacont inui ties. Let repcssent differentiation with reapect to x*
dfx) - D(p(x>) - 0(a(x)) + D(p(x))
Now, D(s(x)) Hill be finite everywhere, while 0(r(x)) Mill be zero aalde
from a number of pulses - which carry ail the information. Each pulae
Mill correspond to an edge between regions and have area proportional to
the intensity atep. If now one ■thresholds'* and discards all finite
parte, one getei
Km) - T(D(p(x)>) ■ D(r(»))
To obtain r(x) one only has to invert the differentiation, that is,
integrate. Let I represent integration with respect to x f then (I) -
D andi
IlKl - I(T(D(p(x))) - rtx) + c
One can give a convolutions! Interpretation to the above, since
differentiation corresponds to convolution with a pulse-pair, one
ON LIGHTNESS 12
negative and one positive, each of unit area. Integration corresponds
to convolution with the unit step function.
1.4.2 Normalization*
The reeult is not unique because of the constant introduced by the
integration. The ztro (spatial) frequency term has been loet in the
differentiation, so cannot be reconstructed. Thle le related to the
fact that one does not knoM the overall level of Illumination and hence
cannot tell whether en object appears dark because it is grey or because
the level of illumination is loti.
One can normalize the result If one assumes that there are no light
•oyrces in the field of view and no floureecent colore or specular
reflections, Thle is certainly the case for the flondrians, Perhape
the beet way of normalising the result Is to simply assume that the
highest value of lightness corresponds to white* or total reflectance in
the Lambertian sense. This normalization will lead one astray if the
image doee not contain a region corresponding to a white patch in the
ecene, but this is the beet one can do. Other normalization techniques
might involve adjusting weighted local averages, but thle would then no
longer amount to reconstruction of reflectance.
ON LIGHTNESS 13
l.S Qne-Oimenaional Hethod - Discrete Case.
So far ue have assumed that tha image intensity uan a continoue
function. In retinas found In animal a or artificial onss constructed
out of dlacrata componsnte. iniagaa ara only saitplsd at discrete points'.
So ona haa to find diacrats analogues for the operations we have been
using. Perhaps tha simplest ar* first differences and summation as
analogues of differentiation and integration respectively. This is not
to say that other approxi mat ions could not be used equally wall.
To use the new operators one goes through essentially tha sane
process as before, except that now all values in the differenced imaga
ara finite. This has the effect of forcing one to choose a threshold
for the thresholding function. Both components of image intensity
produce finite values after tha differencing operation. The component
due to the edges in the reflectance is hopefully quits large compared to
that due to illumination gradient. One has to find a level that Mill
suppress the illumination gradient inside regions, uhi la permitting the
effecte due to edges to reaain,
1*8.1 One-Dimensional Discrete flethod: Details.
Let rl be the reflectance of the object at the point corresponding
to the image point 1. Let s\ be the incident light intensity at this
object point. Let p\ be their product* that la tha lntenelty In tha
image at point I.
- .*
D
j=l
ii
«;
FIGURE 2: Processing steps In the one-dimensional discrete case,
ON LIGHTNESS IS
Now 1st pj b* th* log of p- and bo on. Lot Oft and Ift 0* !"• operators
corresponding to taking first differences and summation respect Ively.
Note that (Ift)"' - D*.
Pi- •l+n
d i- p*h- p; (d - d * ( p ,j
t ; - d; If !d,[ < «. tlio
-£'
ll- > t K -
K-0
1*5.2 Selecting the Threehotd.
Uhat determines the threshold? It mutt be smaller than the
emalleet inteneitg etep between regions. It ituet on the other hand be
larger than values produced by first differencing the maximum
Illumination gradients. Real images are noisy and the threshold ehould
be large enough to eliminate this noise inside regions*
The spacing of the seneor cells must also be taken Into account*
GN LIGHTNESS 16
As this spacing becossv smaller, the contribution due to Illumination
gradients decreases, uhlls the coaponsnt due to the edges remains
constant. A limit is reached uhen the component due to illumination
gradients falls below that due to noise or when the optical properties
of the Imaging system begin to have a deliterlous effect. In all
imaging systems an edge Is spread over a finite distance due to
diffraction and uncorrected abberatione. The spacing of seneore should
not be much smaller than this distance to avoid reducing the component
due to edges In the differenced iaags.
Let u be the radius of the point-spread-function of the optical
system and h the spacing of ths ssnsor celle. Let g* be the smallest
atop in the logarithm of raflectanca In tha acana. Than define tha
affective minimus step asi
g - g* a mind, h/2u)
Let a be the largest slope due to illumination gradient and * the root-
mean-equare noise-amplitude. The noise will exceed a value 3 only «3%
of the time. Choose the threshold e as follows*
a < g
• > a * h
e > 3\/2 a
ON LIGHTNESS 17
1.6.3 Accuracy of the Reconstruction.
In thA contlnoue case one can exactly reconstruct tha reflectance,
aeide from a constant. Us are not eo fortunate here, even If ue Mlect
a threshold according to the above criteria. This la because the valuee
at the edges contain email contributions due to Illumination gradient
and nolee. A alight Inaccuracy In the reconetructlon ulll reeult. Thle
error le minimized by making the esneor cell spacing very fine,
optimally of a elze comaeneurate with the optical reeolutlon of the
device. The effect of nolee can also be minimized by integrating over
tlM.
Note that the reconstruction le aore accurate when there are fen
edgee, elnce It le at the edges that the error effects appear. Ulth
many edgee the Illumination gradient begins to "show through - .
1*5.4 General izat lone.
So far we have dealt ulth constant sensor epaclng. Clearly aa long
ae the aame epaclng le ueed for both the differencing and the summing,
the cell epaclng can be arbitrary and has little effect on the
reconetruction elnce It does not enter Into the equatlone.
■
Similarly ue have chosen first differences as the dlecrete analogue
for differentiation. Ue could have chosen some other uelghted
difference and developed • eultable Inverse for It. Thle Inveree of
course uould no longer be eummatlon but cin be readily obtained using
ON LIGHTNESS 18
techniques developed for dealing with difference equetlone (Rlchtmeyer
1957, Garabedlen 19G4i.
1*5.5 Physical rtodels of the Ons-Oimensional Discrete Proceee.
One can invent a number of physical modele of the above operatlone.
A simple reeistive netuork ulll do for the summation proceee for
example. Land has implemented a small circular "retina" ulth about 1G
eensora. This model employs electronic componente to perform the
operations of taking logarithms, differencing, thresholding and summing.
Land hee tried to extend his one-dimensional method to Images* by
covering the image ulth paths produced by a random ualk procedure and
applying methods like the above to each of these paths. Uhllo this
produces reeulte. It esens unsatisfactory from the point of view of
euggeetlng possible neuro-physlologlcal structureet neither does It
lend itself to efficient Implementation.
Methods depending on non-linear processing of the gradient along
paths in the image fall to smoothly generalize to two dimensions, and
cannot predict the appsarance of images in uhich different pethe reeuit
In different lightneseee.
ON LIGHTNESS 19
2.' Lightness in Two Dimensional [mages.
2.1 Tuo-Oi mansions! flathod - Continous CasSt
Ue naed to extend our Ideas to tuo dimensions In order to deal wl th
actual images, Thsra are a number of Mays of arriving at tha process to
be described here, we shall follow tha simplest (Horn 19681. Urn need to
find two-dimensional analogues to differentiation and Integration. The
first partial derlvatlvee %r% directional and thus unsuitable alnce they
Mill for example completely eliminate evidence of edges running In a
direction parallel to their direction of differentiation. Exploring the
partial derivatives and their linear combinations one finds that the
Laplacian operator la tha louest order combination that la leotroplc, or
rotational ly symmetric. Tha Laplacian operator is of course the sum of
the second partial derivatives.
2.1.1 Applying the Laplacian to a flondrian.
Before investigating the invert ibi I i ty of thie operator* let ue eee
yhat happene when one appliee It to the image of a Mondrlan. Inelde any
region one Mill obtain a finite value due to the variation In
Illumination Intenelty. At each edge one Hi 1 1 get a pulae pair, one
positive and one negative. The area of each pulse ulll be equal to the
Intenelty step.
. -*
<<
1
>
2 j
FIGURE 3: Applying the laplacian operator to the image of a
Hondrian figure.
ON LIGHTNESS 21
Thia can beat be esen by considsring tha first derivative of a
step, nanely a single pulee. If thie ia differentiated again one
obtains a doubled pulse as described. Since this pulse Mill extend
along the edge, one mag think of it as a pulse-wail. So each edge
separating regions m1 1 1 produce a doubled pulse nail. It la clear that
one can once again separata the coaponant due to reflectance and
i I I umi nation simply by discarding all finite parts.
2.1*2 Inverse of the Laplacian Opsrator.
To co«p)ete the task at hand one than has to find a process for
undoing the effect of applying the Laplacian. Again there are a number
of approaches to this problem, we will use the ahortest IHorn 1968). In
essence one has to solve for p(x,y) in a partial differential equation
of the formi
L(p(x,y)) - d(x,y)
Thie ia Poiason'e equation and It la usually solved inside a bounded
region using Green's function (Garabsdien 13541 1
«K,y) - MC{L^« t g)td^^)d{di,
ON LIGHTNESS 22
The form of Gnm'i function G» dependa on the ehape of the region
boundary Nom If the retina ie infinite all polnta are treated
similarly and Green* a function dapenda only on tuo para we tare, (J; - x)
and (*) - y) . Thla poeit tonal indepandanca lapllee that the above
integral almply beconea a convolution. It can ba ehown that Green'e
function for thi a caaa lai
G(£ t *J * *,y) - Q/2 * ) log f (l/r)
Uhere r x - <£,- x) + (*|- y)
So p(x
,y) - / / (1/2") log ( (l/r) * dt^,*| ) d^ d^
Thua the Invaraa of the Laplacian operatora la almply convolution ulth
(1/2*) log <l/r). To ba preciee one haet
(|-=+ i*£ // d/2"> log c (l/r) * d(^,T ) d^dij - d(x,y)
Thie ie the two-dimeneional analogue oft
ON LIGHTNESS 23
*t
fit) dt - f(x)
2*1.3 Why one can use the Convolutional Inverse.
If the retina le coneidsred infinite one can express the Inverse ae
a elmple convolution. If the retina le finite on the ether hand one hae
to uee the acre complicated Creen*8 function foraulation.
Now conelder e scene on a uniform background whose Image le totally
contained on the retina. The result of applying the forward transform
and thresholding will be zero In the area of the uniform background.
The convolutional Inverse will therefore receive no contribution from
outside the retina. As a result one can use the convolutional form of
the Inverse provided the Image of the ecene le totally contained within
the ratine.
2.1.4 Normal Ization.
Once again one flnde that the reconstructed reflectance le not
unlqua. That le. any non-singular solution of Llpfa.y)) - 8 can be
added to tha Input without affecting the result. On the Infinite plane
auch eolutions have the form p(x t y) » (a*x + b*y + c). If the scene
only occupies a finite region of space it can be further shown that the
eolution will be unique up to a constant and that one doee not have to
ON LIGHTNESS 24
worry about possible slopes. To be specifics the background around the
ecene Ml I I be constant in the reconstruction. So one has here exactly
the same normalization problem as In the ons-dleenslonal case.
Assigning uhlts to the region with highest numerical value In the
reconstructed output appears to be e reasonable method.
2*1.6 Tuo-Dimeneional Continous (lethodi Details*
Let r'(x.y) be the reflectance of the object at the point
corresponding to the iaags point <x t y). Let e* (x,y) be the source
intensity at that object point. Let pMx.y) be their product, that ie
the intensity at tht iaage point U.y). Note that rMx.y) and a*(x,y)
are poei tlve.
p* (x,y) - a* (x t y) ft p 1 (x,y)
Let p(x a y) be the logarithm of p* (x.y) and so on.
p(*.y> - a(x,y) t r(x«y)
Now assume that s(x,y) and Its first partial derivatives are continoue -
a reasonable assumption to make for the distribution of illumination on
the object. Let L be the Laplacian opsrator*
ON LIGHTNESS 25
d(x,y) - L(p(x,y)) - Llsfx.y)) + L(r(x,y))
Nom L(i(k»u)I Mill be flnl te everywhere, whl le L(r (x,y) ) mIII bo zero
except at each edge eeparating regions, where one nil I flr>d a double
pulse wall ae deecribed. Now discard all finite partet
t<x,y) - T(L(p<x.y))> -Ltrlx.yM
Lot G bo the operator corresponding to convolution by (1/2* ) log,(l/r]
Note that (G)" 1 - L.
■U,y) - G(T(L(p(* t y))) - r(x t y) +C
2.2 Tuo-Dlmenslonel flethod - Oiecrete Case.
Once again we turn fron a contlnous image to one sampled at
dlecrata points. First Me wit l have to decide on a teaselatlon of the
Image plane*
2.2.1 Teaselatlon of the .nage Plans.
For regular teseelations the choice la between triangular, aquara
and hexagonal unit eel la* In much past work on image proceeding, aquar
teaeelations have been used for the obvious reaeone. Thle particular
taeaalation of the Image has a number of disadvantages. Each cell hae
ON LIGHTNESS 26
two kindt* of neighbors, four adjoining the sides, four on the corners.
This results In a number of aaymmstr lee. It makes It difficult to find
convenient difference schemes approximating the Laplacian operator with
low error term.
Triangular unit cells are evsn worse In that they have three kinde
of neighbors* compounded these drawbacks. Note also that near-circular
objecte pack tightsst in a pattern with hexagonal cella. For these
reasons hs mIM use a hexagonal unit cell. It should be kept In Bind
however that It is easy to develop equivalent results using different
tesselat Ions.
2.2.2 Discrats Analogue of tha Laplacian.
Having decidsd on ths teaselation we nssd now to find a discrete
analogus of the Laplacian operator. Convolution Mlth a central positive
value and a rotational ly syuetric negative surround of equal weight le
one possibility. Aside from a negative scale factor, this Mill approach
application of the Laplacian in the Hull as ths cell size tends to
zero*
If one were to use complicated surrounds, the trade-offs between
accuracy and resolution would suggest using a negative surround that
decreasss rapidly outward. For the sake of simplicity we will choose
convolution with a central cell of weight 1, surrounded by elx eel le of
weight -1/6. Thle function is convenient, symmstrlc end has a small
error ten* It Is equsl to - ( h/4 L + h/64 L 1 ) plus sixth and higher
FIGURE 4a: A discrete analogue of
the laplacian operator.
. -*■
Nflwvtti Delta function minus
this discrete analogue.
ON LIGHTNESS 28
order derivative* IRlchtmeyer 19571 . It should again be pointed out
that similar results can be developed for different functions.
2-2.3 Inverse of the Oiscrete Opsrator,
The forward differencing operator has the form
«»i ■ p., -J}^h * **\
Uhere p,« 1* the logarithm of linage intensity, m ( ; fti'e Might*, uhleh In
our case are 1/G, and the sun is taken over the six immediate neighbors*
Ue now have to determine the inverse operation that recovers d--
from d|i* One approach ie to try and solve the difference equation of
the form
v5j Vi ^ p w " -| i
Or In matrix form U p • d. Note that U la tparse, hevlng l*e on the
diagonal and -1/G'e scattered around. For a finite retina with n eeneor
cells one has to introduce boundary conditions to ensure that one has as
many equations ae there are unknowns. One then eimply inverte the
ON LIGHTNESS 29
matrix U and getst p - U~ d.
This is entirely analogous to the solution in the continous case
for a finite retina. U~ corresponds to the Green's function. Much as
Green's function has a large "support", that is, is non-zero over a
large area, so 1*1 is not sparse. This implies that a Jot of
computation is needed to perform the inverse operation.
2.2.4 Computational Effort and Simplification.
Solving the difference equations for a given image by simple Gauss-
Jordan elimination requires of the order of n /2 arithmetic operations.
Another approach is to invert U once and for all for a given retina.
For each image then one needs only about n x arithmetic operations. Note
that the other operations, such as forward differencing, require only
about Gftn ar i thmet ic operations.
Uhat in effect is happening is that each point in the output
depends on each point in the differenced image. Both have n points, so
n 1 operations are involved. Not only does one have to do a lot of
computation, but must also store up the matrix U of size n . This le
quite prohibitive for even a small retina.
This latter problem can be avoided if one remembers the
simplification attendant to the use of an infinite retina in the
continous case. There ne found that the integral ulth Green's function
simplified into a convolution. Similarly, if one assumes an infinite
retina here, one finds that W and its inverse become very regular, Tho
ON LIGHTNESS 38
rows in U are then all the sane and the same is true of U . Each value
in the output then depends in the same way on the neighboring points In
the differenced image. One need only store up the dependence of one
point on its neighbors for this simple convolutional operation.
The only remaining difficulty is that U 19 now Infinite and one can
no longer invert it numerically - one has to find an analytical
expression for the inverse. I have not been able to find this inverse
exactly. A good first approximation is log,lr D /rl - except for r » 0,
when one uses 1 + log,(r c K Here r is the distance from the origin and
p^ is arbitrary. The reaainder left over i-ihen one applies the forward
difference scheme to this approximation lies between log.ll + r ) and
log 11 - p" ) t This error term is of the order of r •
In practice one does not have an infinite retina, but as has been
explained for the continous case one can use the convolutional method
described above for a finite retina, provided that the image of the
scene is wholly contained inside the boundaries of the retina* It is
possible to find an accurate inverse of this kind valid for a limited
retinal size by numerical means.
2.2.5 Two-Dimensional Discrete Method: Details.
Let r*. be the reflectance at the object point corresponding to the
image point (l.j). Let »y be the intensity of the incident light at
this object point. Let &l be the intensity in the image at point
ON LIGHTNESS 31
»?j ■•b-4
Let p;j be the logarithm of pV and so on. Ltt L* be the operator that
corresponds to convolution with tha analogua of tha Laplaclan* Lat G*
be i ta Inverse.
df i * p 'i "2-?^* p w ,d - u( p» J
Tha weigthe m- are 1/6 in thte casa t and tha eua la taken over the ah
Immediate nelghbora.
else
Here the sum extends over the whole retina and v-» is the convolut ional
Inveree found numerically as explalnsd above.
ON LIGHTNESS 32
2.2.G Simplicity of the Invert*.
The forward transform. Involving only a slmpls subtraction of
Immediate neighbors, Is clearly a rapid, local operation. The inverse
on the other hand ie global, since each point in the output depends on
each point in the differenced image. Computationally this makes the
inverse slow. The inveree is simple in one sense houeveri The
difference equations bslng solvsd by ths Inverse have the seme form as
the equations used for the forward transform end are thus local. The
problem Ie that the output here feeds back into the system and effecte
can propagate across ths retina. The apparent global nature of the
inverse is thus of d rather special kind find, as we will see later,
gives rise to very simple implementations involving only local
connect tone.
2.2.7 Iterative Methods of Solution.
There are of course other methods for solving large sete of
equations. The fact that U is sparse and has large diagonal elements*
euggeets trying something I ike Gauss-Ssldsl iteretion. Each Iteration
takes about S*n arithmetic operations. For effects to propagate acroee
the retina one requires at least ^(4ftn* - l)/3 iteratlone. Thle le
becauee a hexagonal retina of width m has (3*m x + l)/4 cells. The above
euggeets that one might be able to get away with less than n l arithmetic
operatione. In practice It Is found that effecte propagate very elouly
ON LIGHTNESS 33
and many more Iterations art needed to stabilize the solution, Ona does
not hava to atora U t alnca it la aaailg generated aa one goea along.
Iterative echemee correspond to adding a tims-der ivative to tha
Poiaaon equation and so turning It into a heat-equation, Aa ona
continues to iterate the steady-state solution la approached. This
intuitive modal givea some insight into how the process will converge*
2.2*8 Convergence of Iterative and Feed-back Schemes.
It is not immediately clear that iterative achamea of solving the
difference equations hIJI converge. If they do, they will converge to
the correct solution. Let 5 be the delta function, that is. one at tha
origin, zero elsewhere. It can be shown that if tha forward
convolutional operator la u* the convergence of Iterative echemee
depends on the behaviour of the error term, ( 6 * w) , aa n bacomaa
large. Raiaing a convolutional operator to an integer power is
intended to elgnlfy convolution with itself.
In our case, m is ons at the origin, ulth elx valuea of -1/6 around
it. So ( 6 - y) ill 1 1 b« zero at the origin with elx valuaa of 1/6 around
it. New while ( i- w) Ml 1 1 aluays have a total area of one. It doea
epread out and I te value tends to zero at every point aa n tende to
Infinity. So thie iterative scheme convsrgeet similar results could bo
derived for other negative surrounds.
ON LIGHTNESS 34
2.2*9 Setting the Threshold.
In the dlecrete case a finite threshold must be selected. As
before, let g' be the seal lest stsp In the logarithm of reflectance In
the ecene, h the sensor spacing and u the radius of the point-spread
function of the optical system* Than us define the effective minimum
step asi
g • g' * mind, h/2u)
There are some minor differences In what follows depending on whether
one considers the sensor outputs to bo intensity aanplee at cell-centore
or averages over the cell area. The smallest output dus to Mn edge ui II
be about g/6. This Is producsd uhsn ths edge le oriented to cover Just
one cell of the neighborhood of six. Let Q be the maximum of the
intensity gradient - that is the Laplacian of intensity In this case.
Choose the threshold s as follows:
e < g/G
e > * h*
e > 3^7/G" *
ON UGHTNESS 35
2*2*10 Sons Notes on This tlethod*
Notics that an illumination gradient that varies aa soms power of
distance across ths image becomes a linear slops after taking logarithms
and thus produces no component after the differencing operations* Such
simple gradients are eupprseeed even without the thresholding operation.
In practice ths parameters used in choosing the threshold may not
be known or may be variable. In this cass ons can look at a histogram
of the differenced image. It will contain values both positive and
negative corresponding to edges and also a large number of values
clustered around zsro due to i I lunination gradients, noise and so on*
The threshold can bs convsnisntlu choosen to contain this central blob*
Noise and illumination gradients havs an effect similar to that In
the one-dimensional cass. Ulth finite cell spacing ons cannot precisely
eeparate ths two components of ths image intensity and at each edge the
information Mill be corrupted slightly by nolss and Illumination
gradient. As the density of sdgas per csll arsa goss up the effect of
this becomes mors apparent. In highly textured scenes the illumination
gradient is hard to eliminate*
Once again one hee to dsclde on a normalization echeme. The best
method probably Is to 1st the highest numerical value In the
■
reconstructed output correspond to unite*
ON LIGHTNESS 36
2.2.11 Dynamic Range Reduction.
Applying the retinex operation to an image coneiderable reducee the
range of vatuee. Th'io ie because the output, being related to
reflectance, will only have a range of one to two ordere of magnitude,
while the Input will a I eo have i lluminat ion gradients. This will make
euch processing useful for picture recording and transmission (Horn
19681*
2.2.12 A Frequency Domain Interpretation*
t
It may be of Intoreftt to look at thit as t hod fron got another point
of view. Uhat one does Is to exentuate the high-frequency components,
threshold and then attenuate the high-frequency coeponente. To eee
this, coneider first the forward operation. The fourier traneform of
the convolut tonal operator corresponding to differentiation le 1«.
Similarly the two-dimensional fourier traneform of the convolutional
operator corresponding to the LapJeclan le - p*. Here p ie the radiue
in a polar coordinate system of the two-dimensional frequency space. In
either case one is multiplying the fourier transform by sone function
th^t increaeee with frequency. Now consider the revsree operation. The
fourier traneform of the convolutional operation corresponding to
Integration ie l/l«. Similarly the fourier transform of (1/2*) log r
(1/r) ie -l/p x . So In the Inverse step one undoes exactly the emphaele
given to high frequency components in the forward operation.
ON LIGHTNESS 37
In both thn ons-dlmenslonal and the tuo-dinensional case one loeee
the zero frequency component. Thle It why the reeult hae to be
normal ized.
2.3 Phgelcal llodele.
There are nueerous continous physical nodele to illustrate the
inverse traneformatlon. Anything that satisf lee Polsson'a equation ulll
do. Such physical models help one visualize Mhat the Inverse of a given
function might be. Examples in two dimensions arei perfect fluld-f Iom»
eteady diffusion, steady heat-flou, deformation of an elastic membrane,
electro-statics and current flow in a resistive sheet. In the last
model for example, the input is the distribution of current flowing Into
the reeietlve sheet noraal to Its surface, the output ie the
dletrlbutlon of electrical potential over the surface.
In addition to helping one visualize solutions* these contlnoue
eodels aleo suggest discrete models. These can be arrived at elaply by
cutting up the two-dimensional space in a pattern corresponding to the
interconnection of neighboring cells. That Is, the remaining parte form
a pattern dual to that of the sensor cell pettern. Us ulll discuss only
one such discrete mods It
ON LIGHTNESS 38
2.3.1 A Discrete Physical node I.
Consider the reeletive eheet described* cut up in the dual pattern
of the hexagonal unit cell pattern. What will be left ie an
Interconnection of resistors In a triangular pattern. The inpute to
thie system will be currents injected at the nodes, the potential at the
nodes being the output. Thie then provides a very simple analog
implementation of the tedious inverse computation.
It ie perhaps at first surprising to sss that each cell ie not
connected to every other In a direct fashion. One would expect thie
from the form of the computational Inverse. Each cell In the output
does Of COUrft* h«V* ft connection Via the ftthtr ceils to oich of thi
Inputs. Paths are shared however in a way that makes the result both
simple and planar.
Conelder for the moment just one node. The potential at the node
ie the average of the potential of the six nodes connected to it plue
the current injected times H/B, where R Is the resistance of each
resistor. The economy of connection is due to the fact that the
outpute of thie system are fed back into it. It also illuetrates that
thie model locally solves exactly the same difference equation ae that
used In the forward transform, only now in reverse.
This immediately suggests an important property of thie model! By
eimply changing the Interconnections one can make an inverse for other
forward transforms. Slmpleet of all are other image plane tesse lat ions,
both regular and irregular. One simply connecte the resistors in the
■ '*
FIGURE 5: Resistive model of the inverse computation.
The inputs are the currents injected at the nodes.
The outputs are the potentials at the nodes*
ON LIGHTNESS 40
i
sane pattern as are tha calls in tha input*
Mora complicated weighted surrounds can ba handled bg using
resistors with resistances inversely proportional to tha weights* Tha
network of resistors will then no longer be planar.
2*3.2 A Feed-back Scheme for the Inveree.
Both the comment about outputs feeding back into the resistive
Model and the earlier notes about iterative schemes suggest get another
interesting model for the inverss using linear cunning devices.
Operational amplifiers can serve this purpose. One simply connects the
summing element so that they solve the difference equation implied by
the forward transform. Once again It la clear that such a scheme can be
generalized to arbitrary tesselations and wsighted negative eurrounde
eimply by changing the Interconnections and attenuations on each input.
Some questions of stability arise with esotsrlc Interconnections. For
the simple ones stability Is assured.
A little thought will show that ths resistive model described
earlier ie in fact a more economical implementation of just this scheme
with the difference that there the inputs are currents, while here they
are potentials.
£
u
*j
id
.a
TJ
c
i
m
"E
■u
2
c
t.
<0
U1
4-t
e
2
i
+J
a
3
c^
■
c:
M-
*C
o
c
£
^
o
VI
40
•*—
c
*>
<v
"-
^
*-
o
*-*
ft)
8
V
^
^
^
-j
1-
CJ
I
>
*^
ON LIGHTNESS 42
2.4 Limitations of the Simple Scheme Presented.
Tha Method presented hara ul II not correctly calculate reflectance
if uaed unmodified on general scenes. It may however calculate
lightness fairly well. Ae the msthod stands now for example, a sharp
ehadow edge ulll not be distinguished from a real edge In the scene and
the two regions so formed nil 1 produce different outputs, whi le their
reflectances are the same. It nay be that thie le reasonable
nevertheleee, elnce ue perceive a difference in apparent lightness.
Smooth gradations of reflectance on a surface due either to shading
or variations In surface reflectance ulll be eliminated by the
thresholding operations except as far as they affect the intensity at
tha borders of the region* Thie may imply that ue need additional
channele In our visual eystem to complemsnt the onee cerrylng the
retlnened Information elnce ue do utilize shading as a depth-cue*
The simple normalization scheae described Mill also be sensitive to
specular reflections, flourescent paints and light-sources in the field
of view* Large depth-discontinul ties present another problem. One
cannot assume that the illumination is equel on both sides of the
obscuring edge. In this case tha illuminating component does not vary
smoothly over the retina, having inetead some sharp edges.
S -2
C ** *->
3 £ E
as
.woe
E ^
0> r-
5 2
'2 i
E 5
4
5
V
a
>
in
o
Id
CI
a.
o
L/l
4B
3
It
1%
o
L3
ON LIGHTNESS 44
2.S Computer Simulation of the Discrete Method*
A computer program was used to simulate the retinex process
described on a small retina with both artificial and real images eeen
through an image dissector camera. The hexagonal unit cell le used in
this program and the retina itself Is also hexagonal. The retina
containe 1027 cells in a pattern 36 celle across. Thie ie a compromise
dictated by the need to limit the number of arithmetic operatlone In tha
inveree traneform. In this case one needs about a million and thie
takee about a minute of central processor time on our PDP-18.
Both the art 1 f Iclal and the real Hondrians consist of regions
bounded by curved outline* to smphaaiza that this ntthod does not
require straight-line edges or boundsry extraction and description.
Various distributions of Incident illumination can be selected for the
artificial ecenes. In each case the processing satisfactorily removes
the gradient.
For the real scenes It is hard to produce really large Illumination
gradiente by poeltloning the light-sources. The reconstruction does
eliminate the gradient well, but often minor flaws m1 1 1 appear In the
output due to noiee in the input and a number of probleme ulth this kind
of Input device such as a very considerable scatter. It ie not easy to
predict uhat effects euch imaging device defects will have.
The output is displayed on a OEC 349 display which has a mere eight
grey-levele. It uould be interesting to experiment ulth larger retinae
and better Image Input- and output-devices.
ON LIGHTNESS AS
2.5.1 For« of Inverse used In the Computer Simulation*
The convo I u t i one I form of the Inverse was used for speed and lou
storage requirement. This neccesitated solving the difference equation!
once, given a pulse as Input. The symmetry of the hexagonal pattern
allous one to identify symmetrically placed cells and only 324 unknowns
needed to be found for a convolutional inverse sufficient for the eize
of rotina described. As Mentioned before, this function is closely
approximated by log, (r # /rl for large r. This can be used to establish
boundary condi t ions.
FIGURE 8: The method applied to an
artificial image.
FIGURE 9: The method applied to a
real image
FIGURE 10: The method applied to
Craik's figure.
FIGURE 11: Apparent lightness predicted
for incomplete figure-
The subfigures in the above have the following interpretation:
A Input - logarithm of image intensity
B Differenced image
C Thresholded difference
D Output - computed lightness
E Illumination distribution (P^*!**)
1
fj
■i.i
ON LIGHTNESS 47
3* Implications and Conclusions.
3*1 Parallel Imaga Procsaaing Hardware.
Tha methods described hare for forward transforming, thresholding
and Inverse transforming inaediately tempt one to think In terms of
electronic componenta arranged in parallel layere. Enough has been aald
about different models to make it clear how one Bight connect euch
components. Large scale integrated circuit technology may be useful*
provided the signals are either convsrtsd from analog to digital for* or
batter still, good linear circuits are available In thia form.
Construction of such devlcea would be premature until further
experimentation le performed to decide on optimal tesselatlone, optimal
negative surrounds, thresholding operations and normalization schemes.
Theee decleions are best guided bg computer simulation.
3.2 Cognitive Peychology.
One of the artificial ecanes was created to Illustrate Cralk's
lllueion (Brindley I960, Cornsweet 1970). Here a sharp edge le bordered
by eecond-order gradients. As one might expect the smooth gradlente are
loet in the thresholding and reconstruction produces two regions each of
uniform brightness* Ths diffsrencs in brightnsss between the reglone le
equal to the original Intensity stsp at the edge.
The fact that the process presented hsre felle prey to this
ON LIGHTNESS 48
i I lueion is of course no proof that humans use the «a»* mechanism. It
le interesting that this technique allows one to predict for example the
appearance of pictures containing incompletely closed curves with
second-order gradients on either side.
3*3 Neuro-physr ology.
The method described here for obtaining lightness from Image
Intmnelty euggeete functions for a number of etructuree in the primate
retina. The horizontal ceils appear to be involved In the forward
transformation, while some of the amacrlnes may be Involved in the
inverse transformation* For details see tha paper by David Narr IDarr
19741 , In which he uses this hypothssis to explain an astonishing number
of facte about the retina.
3.4 Concluelon.
A elmple layered, parallel technique for computing lightness from
image Inteneity has been preesntsd. Ths method does not involve an
ability to describe or understand the ecene, relying instead on the
spatial differences In the distribution of reflectance and illumination.
The forward step involves accentuating ths edges between regions. The
output of thle step is then thresholded to remove Illumination gradiente
and noise. The Inverse step merely undoes the accentuation of the
edges*
ON LIGHTNESS 49
Physical models nave been glvan uhich can perform this commutation
efficiently in parallel layers of simple networks. The method haa been
simulated and applied to a number of images* The method grew out of an
attempt to extend Land'e method to two dimenslone and fills the need for
a llghtneee Judging process in his retinex theory of color perception.
The possibility of processing an imags In such a parallel* elmpla
fashion uithout higher-level understanding of the scene reinforces my
belief that such lou-level processing is of Importance In dealing ulth a
number of features of Images. Amongst thsee are shading, stereo
disparity, focue, edge detection, scene segmentation and motion
parallax. Some of thia kind of processing mag actually happen In the
primate retina and viaual cortex. The Implications for Image analysis
are that it may well be that a number of such pre-proceeeing operatlona
should be performed automatically for the whols Image to accentuate or
extract certain attributes before one brings to bear the more powerful*
but tedious and slou sequential goal-directed methods.
CN LIGHTNESS SB
B lb I iographg.
Brindleu. G.S, (1368) "Physiology of the Retina and Visual Pathway."
(Honograph No. 6 of the Physiological Society). Londoni Edward
Arnold Ltd.
Cornsueet, T. (1978) "Visual Perception." Nau Vorki Academic Preae.
Garabedian, O.R. (19G4) 'Partial Olffarantlal Equations." Nau Yorki John
Ui lay.
Goetha, J. 14. von. (1818) "Zur Farbenlehre. " Tuebingen.
Hardy, A.C. (Ed,) (193S) "Tha Handbook of Coioriaetry. " Cambridge, Haest
n. I ,T» Praaa*
Hefmhoitz, H.L.F. (1867) "Handbuch dar Physlologlachan Optlk." Lalpzlgi
Voaa. Alao tranalatedi Southall, J.P.C.S. "Handbook of
Physiological Optics." Nau Vorki Oovar Publications.
He I son, H. (1938) "Fundamental Problems In Color Vision I." Journal of
Experimental Psychology, 23.
Helaon, H. (1948) "Fundamental Problems in Color Vision II." Journal of
Experimental Psychology, 26.
Herlng, E. (1875) "Zur Lahrs von Lichtsinne - Grunzuege elner Theorie
des Farbsinnes." SBK. Akad. Uiss. Uien ftath. Naturuies. K.7B.
Horn, B.K.P, (1368) "Tha Application of Fourier Transform Methods to
Image Procasalng." S.M. Thesis. E.E. Department, fl.I.T. pp. 81-87 A
pp. 93-94.
Judd, D.B* (1948) "Hue, Saturation and Lightness of Surface Colore with
Chromatic Illumination ." Journal of the Optical Society of
America, 38.
Judd, D.B. (1962) "Color In Business, Science and Industry." New Vorki
John Ulley.
Land, E.H. (1959) "Experlmsnte in Color Vieion," Scientific American.
Land, E.H, (1964) "The Retinex," American Scientist. S2.
Land, E.H. 4 flcCann, J.J. (1971) "Lightness Theory." Journal of the
Optical Society, 61.
Lettvin, J.Y. (1967) 'The Color of Colored Things." Quarterly Progress
Report, 87, Rssearch Laboratory for Electronics, M.I.T.
ON LIGHTNESS 51
Marr f D- (1974) "An Analysis of the Primate Retina." A.I. Clemo, 296,
Artificial Intel I Igence Laboratory, M.I.T.
tlaxueH, J.C. (185G) "On the Unequal Sensitivity of the Foramen Centrale
to Light of Different Colours. " Rep. Brit. Assoc.
Newton, Sir Isaac. (1784) "OpticKs." Londoni Sanuel Smith 6 Benjamin
Uatford. Also Neu Yorkt Dover Publ icatione.
Richards, U. (1971) "One-stage flodsl for Color Convereion. ■ Journal of
the Optical Society, 62.
Rlchtmeyer, R.O. & dor ton, K.U. (1957) "Difference Hethode for Initial
Value Problems.- Neu YorKi John U isy.
Young, T. (1828) "On the Theory of Light and Colour." Philosophical
Transact lone* Aleo im Teevan, R.C. & Birney, R.C. (Ed.) (1361)
"Color Vision," Van Noetrand.